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POWER CALCULATIONS FOR ISENTROPIC COMPRESSIONS 


OF CRYOGENIC NITROGEN 

Jerry B. Adcock and Marilyn E. Ogburn 
Langley Research Center 


SUMMARY 

A theoretical analysis has been made of the power required for isentropic 
compressions of cryogenic nitrogen. This real-gas analysis was made from a 
cryogenic wind-tunnel perspective, and the purpose of the analysis was to deter- 
mine the extent to which wind-tunnel drive power might be affected by the real- 

gas characteristics of nitrogen. The real-gas solutions covered stagnation tem- 
peratures from 80 to 310 K, stagnation pressures from 1.0 to 8.8 atm, and fan 

pressure ratios from 1.025 to 1.200. These solutions are compared to the ideal 

diatomic gas solutions. At cryogenic temperatures, the power required to com- 
press nitrogen isentropically is as much as 9.5 percent lower than that required 
for the ideal gas. Simple corrections to the ideal values of mass flow, energy, 
and power were found to give accurate estimates of the real-gas values. 


INTRODUCTION 

The cryogenic wind-tunnel concept has been developed at the Langley 
Research Center in order to improve flight simulation in wind tunnels by 
increasing the test Reynolds number. The major advantages of increasing the 
Reynolds number by reducing the temperature of the test gas are given in refer- 
ences 1 to 4. For fan-driven tunnels, one of the prime advantages of reducing 
the temperature of the test gas is the accompanying reduction in the required 
drive power. The extent of this reduction is illustrated in figure 1 for three 
different wind-tunnel cases. In each case two of the three tunnel parameters 
(Reynolds number, stagnation pressure, and size) are held constant while the 
third varies with decreasing stagnation temperature. These calculations are 
based on the assumption of an ideal gas. However, for the cryogenic wind-tunnel 
concept as developed at Langley, cooling is accomplished with liquid nitrogen; 
the resulting test gas is cryogenic nitrogen. At these conditions, nitrogen 
has real-gas imperfections (ref. 5). Even though the analysis of reference 5 
indicated that cryogenic nitrogen would be an acceptable test gas in terms of 
flow simulation, the real-gas characteristics could possibly become important 
in a wind-tunnel design consideration such as the drive power requirement. 

This report presents the results of a study to determine the extent to 
which the wind-tunnel drive-power requirements might be affected by the real- 
gas characteristics of nitrogen. In this study, real-gas solutions for 
isentropic compressions of nitrogen were made for the range of operating 
temperatures (saturations to 310 K) , pressures (1.0 to 8.8 atm), and tunnel 
pressure ratios (1.025 to 1.200) anticipated for fan-driven transonic cryogenic 



wind tunnels. The' solutions were compared with those for an ideal diatomic 
gas, and the results are presented relative to the ideal-gas values. Real-gas 
calculations for tunnel throat conditions and for fan calculations are given 
in appendixes A and B. Program listing is given in appendix C. 


SYMBOLS 


constants 


speed of sound, m/sec; W in computer program 

specific heat at constant pressure, J/mole-K; CP in computer 
program 

specific heat at constant volume, J/mole-K; CV in computer 
program 

energy per unit mass, J/kgm 


specific enthalpy, J/mole-K; H in computer program 
Mach number 

mass-flow rate per unit area, kgm/sec-m 2 
power per unit area, J/sec-m 2 

pressure, atm (1 atm = 1.013 * 10^ Pa); P in computer program 

gas constant for nitrogen, 296.791 J/kgm-K 

pressure ratio, Pt,2/Pt,1 

entropy, J/mole-K 

temperature , K 

velocity, m/sec; VEL in computer program 
specific volume, liter/mole 
compressibility factor, pv/RT 
specific heat ratio, Cp/c v 
density, mole/liter; D in computer program 
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Superscript: 

a 

a isentropic expansion coefficient, p = P (Constant) 

Subscripts : 

1 upstream of fan 

2 downstream of fan 

amb ambient temperature, 310 K 

i ideal-gas value 

r real-gas or nitrogen value 

S constant entropy 

TH tunnel throat 

TS test section 

t stagnation condition 

u increment 


BASIC EQUATIONS 

The test gas of a closed-circuit fan-driven wind tunnel is forced to flow 
around the circuit by the energy which is imparted to the gas by the fan. If 
the steady-flow compression that takes place at the fan is assumed to be a 
reversible adiabatic process (that is, isentropic), the energy per unit time, 
or power, which must be imparted to the gas is given by the equation 


p = “( h t,2 - h t ,i)s 


( 1 ) 


This equation holds for any gas and thus will be termed the real-gas power 
equation for isentropic compressions. This equation appears to be simple 
enough, but the two factors are not easily calculated. A later section will 
describe the method used in solving this real-gas equation. 

The assumption of an ideal gas and the resulting expressions for isen- 
tropic flow allow this power equation to be expressed in an easily calculated 
form. An ideal gas as defined in this study is one that is both thermally and 
calorically perfect. The ideal-gas characteristics are: 
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1 


pv 

RT 


which is the equation of state, 


J dh = c p JdT 


which gives the specific heats independent of T and p, and 


In addition, the expressions that relate the static variables in an isentropic 
process are 

Y 

p = p Y (Constant) = T Y_1 (Constant) 


By using these characteristics and expressions, equation (1) can be put into 
the following form: 



This ideal-gas power equation and the real-gas power equation ( 1 ) have 
been arranged to show two distinct factors. The first is the mass- flow rate, 
and the second represents the energy per unit mass for each case. This study 
analyzes the real-gas effects on the power required for isentropic compressions 
by examining the manner in which each of these factors is affected. 

The mass- flow rate of any gas per unit area is given by 


m s pV 
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A subsequent section and appendix A describe the real-gas solutions for tunnel 
mass- flow rate. For isentropic flow of an ideal gas, this equation can be 
expressed as a function of the stagnation conditions and Mach number: 


Y+1 



ANALYTICAL MODEL OF TUNNEL 

The real-gas effects on the power required for isentropic compressions of 
nitrogen could be analyzed by assuming that the compressions take place in a 
constant-area duct where nitrogen is flowing at various temperatures, pressures, 
and velocities. However, since the impetus for this study evolved from the con- 
sideration of the power required for the operation of transonic cryogenic wind 
tunnels, the analysis will instead be made from this perspective. 

A sketch of the analytical model of the tunnel is shown in figure 2. For 
this analysis the conditions upstream of the fan p^ i and Tt i and the 
pressure ratio r across the fan are the same for ttie real-gas ’and the ideal- 
gas cases. This assumption means that the outlet temperature T^ 2 is differ- 
ent for each case. With the assumption of no energy losses between the fan 
outlet and the tunnel throat (explained later), the stagnation conditions are 
the same at both of these locations . As a consequence , the test section or 
throat temperature ^t,2 is not the same f° r the real- and ideal-gas cases. 

It will be shown later’ that the difference between (T^^r and (Tt ,2) ± is 
insignificant; therefore, for all practical purposes, the comparisons between 
the real and ideal gases are made at the same test-section conditions. 

The tunnel mass flow will be calculated for the throat conditions (that 
is, Tj. 2» Pt 2» an< * %h). For su bsonic speeds, the throat and test-section 
Mach numbers af*e assumed to be identical. For supersonic speeds, the effective 
area of the test section has to be larger than that of the throat. In practice, 
this larger effective area is created either by diverging the walls of the test 
section or by allowing some of the mass of gas to flow through porous or slotted 
sections of the wall into the plenum chamber. In this latter case, the mass 
may be removed from the plenum by auxiliary suction, or it may reenter the test 
section at the diffuser entrance. For the analytical model used in this study, 
all the mass that passes through the tunnel throat is assumed to pass through 
the fan also. 

For simplification, all the tunnel energy losses are assumed to occur 
between the throat of the tunnel and the fan. Data from existing transonic 
tunnels indicate that most of the losses do occur in this part of the tunnel 
because of the higher flow velocities. However, the results of this report 
could be used for other loss distributions if the fan energy and the throat 
mass flow were evaluated at the appropriate stagnation conditions. 
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For the cryogenic wind-tunnel concept, cooling is accomplished by evap- 
orating liquid nitrogen in the tunnel flow, and as a result, mass is added to 
the stream at the cooler. It is anticipated that the cooling system would be 
placed upstream of the fan rather than downstream because the longer distance 
to the test section would permit more thorough mixing of the evaporating nitro- 
gen with the main stream. This additional mass flow due to cooling is at most 
about 2 percent of the tunnel mass flow. A brief analysis which took this 
additional mass flow into consideration indicated that although the absolute 
level of power was up by 2 percent due to compression of this additional mass, 
the ratio of the real to the ideal power requirement was not significantly 
affected. Thus, for simplicity, cooling of this analytical model is assumed 
to occur without mass addition (mass-flow rate is identical at all tunnel 
locations). 

Typical values of the fan pressure ratio which are necessary to achieve 
a given Mach number in the test section have been assumed for this analytical 
tunnel : 


m ts 

r 

0.2 

1 .025 

.6 

1.050 

1 .0 

1.100 

1.2 

1.200 


This Mach number pressure ratio correspondence is further assumed to. be invar- 
iant with stagnation temperature and pressure. 

This analytical model is , assumed to have an operating stagnation pressure 
range of from 1.0 to 8.8 atm* The maximum pressure matches that of the proposed 
National Transonic Facility that is currently being designed (ref. 6). The stag- 
nation temperatures cover the range from 310 K down to the saturated vapor tem- 
perature. Specifically, this lower limit of stagnation temperature at a given 
stagnation pressure is taken to be that temperature which causes the static 
temperature and pressure at the tunnel throat to be coincident with a point on 
the vapor pressure curve. 


PROCEDURE FOR ANALYTICAL SOLUTIONS 

Figure 3 shows a flow chart of a program that was written in order to cal- 
culate the power required for isentropic compressions of nitrogen and an ideal 
diatomic gas. This program uses a nitrogen properties program written at the 
National Bureau of Standards (ref. 7) that is based on Jacobsen's equation of 
state (ref. 8) It also makes use of some of the subprograms and procedures that 
were developed for the isentropic expansion study of reference 5. Appendix C 
gives a program listing and sample output. 

As the flow chart shows, the program inputs are the test-section stagnation 
pressure Pt,2> Mach number M-pg* and the fan pressure ratio r. First, the 




program sets the test-section stagnation temperature for the real-gas or nitro- 
gen case to a value of 310 K. Next, the program makes two sets of real-gas cal- 
culations. The first of these is the real-gas calculation of the tunnel throat 
conditions (block A). After the throat Mach number has been set, this routine 
determines the static flow properties which would result in the desired M.™. 
When this procedure is completed, the real- gas mass- flow rate is determined. 

As these calculations are being made, a check is made to see whether the static 
flow properties have reached the saturated condition. If saturation occurs, 
the solutions are terminated. The details of the throat calculations of 
block A can be found in appendix A. 

The other set of real-gas calculations is related to the fan as indicated 
by block B. Assuming isentropic compression, this routine takes the downstream 
stagnation conditions p^ 2 and (t^. 2 ) r and the fan pressure ratio and com- 
putes the upstream stagnation conditions p^. ^ and Tj. ^ . The real-gas energy 
per unit mass E r is also determined and combined with* the mass-flow rate from 
block A to give the real-gas power P r for the compression. The details of 
these fan calculations are given in appendix B. 

The next step in the program is to make the fan calculations for the ideal- 
gas case. The inlet stagnation conditions as determined from the real-gas 
calculations and the fan pressure ratio are used in these calculations. The 
outlet temperature is determined from the following isentropic relationship : 


Y-1 


( T t , 2) i = r Y T t>1 


The energy per unit mass is given by 



which is the second factor of basic equation (2). 

After the ideal stagnation temperature (T. 2)^ has been determined, the 
ideal mass- flow rate at the throat is determined 1 from basic equation (3). This 
mass- flow rate is combined with the energy per unit mass to give the ideal 

power required for the compression. 

At this point, the program prints out the real- and ideal-gas parameters 
associated with the compression and their relative values. With the completion 
of this solution, the real-gas temperature T^. 2 is decreased and the solution 
repeated. The temperature is incrementally decreased in this manner until the 
throat conditions for the real-gas case become saturated. 
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ANALYSIS OF SOLUTIONS 


Isentropic power solutions of the type just described have been made for 
the analytical model of a cryogenic tunnel. The solutions cover a test-section 
Mach number range and a fan pressure ratio range from 0.2 to 1.2 and 1.025 
to 1.200, respectively. Solutions covering these ranges were made at stagnation 
pressures from 1.0 to 8.8 atm and at stagnation temperatures from 310 K to sat- 
uration temperatures. This analysis examines the real-gas effects on the power 
for isentropic compressions by showing the effects on the two factors, mass- flow 
rate and energy per unit mass, that combine to give the power. 

It should be remembered that the real- and ideal-gas solutions are for the 
same fan inlet conditions p^ f i and T^ i and for a given pressure ratio r. 

The outlet pressure p*. 2 i s the same, but the outlet temperatures are differ- 
ent. This difference i!k shown in figure 4 for the conditions which produce the 
maximum difference r = 1.20 and Pt 2 = atm. Even at the lowest inlet tem- 
perature, the real-gas value of downstream temperature differs from the ideal- 
gas value by less than 0.3 percent. Thus, for all practical purposes the 
following comparisons of real- and ideal-gas solutions are made at identical 
test-section stagnation conditions as well as for identical fan inlet conditions. 


Energy for Isentropic Compressions 

The real-gas effects of nitrogen on the energy per unit mass for isentropic 
compressions are shown in figure 5. The relative values (real to idrfeal) of 
energy are presented as a function of fan outlet temperature T t 2 anc * ^ or 
various values of outlet pressure p^. 2 . Each curve is for a given fan pressure 
ratio. The fan outlet conditions Pt*2 ' an< 3 Tt 2 were chosen as the independ- 
ent variables because these are the values for the 1 tunnel throat and test 
section. 

These figures show that the nitrogen values for energy per unit mass are 
always less than the ideal values. This difference increases as temperature 
is reduced. At the maximum pressure (fig. 5(d)), the real energy per unit mass 
is as much as 17 percent lower than the ideal diatomic gas value. These lower 
values of energy per unit mass for nitrogen could have been anticipated by 
comparing the curves for enthalpy against temperature at constant entropy 
(fig. 6). For the ideal gas, enthalpy is a function of temperature only. 

Along an isentrope, the enthalpy of nitrogen is a function of both temperature 
and pressure . Since the temperature dependence is dominant , however , the slope 
difference at a given temperature should be an indication of the energy ratio 
for the two cases. The slope for the nitrogen isentrope is less than that for 
the ideal gas. 

Figure 5 also shows that the value of fan pressure ratio has an insignif- 
icant effect on the shape of the energy ratio-outlet temperature curve. 

Figure 7 shows the effect of pressure on the energy ratio at constant 
temperatures. As can be seen, the energy ratio decreases nearly linearly with 
increasing pressure. 
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Real-Gas Effects on Tunnel Mass Flow 

The relative values (real to ideal) of the tunnel mass flow are shown in 
figure 8 as a function of stagnation temperature for various values of stagna- 
tion pressure. As mentioned previously (fig. 4), the outlet stagnation temper- 
atures were so near the same value for the real- and ideal-gas cases that these 
comparisons are at essentially the same throat conditions. Each curve is for 
a given throat Mach . number . 

As stagnation temperature is reduced, the mass-flow rates for nitrogen 
become increasingly greater than those for an ideal diatomic gas. At the maxi- 
mum pressure (fig. 8(d)), the nitrogen mass-flow rate for Mth ="0.20 and the 
minimum temperature is about 9.5 percent greater than the ideal-gas mass-flow 
rate . For Mth = 1 • 0 , the real mass flow is only about 7 . 0 percent greater due 
to the saturation temperature being higher than for the 0.2 case. The shape 
of the mass-flow — temperature curve is relatively insensitive to the throat Mach 
number . 


Power for Isentropic Compressions 

The real-gas effects on the power required for isentropic compressions are 
shown in figure 9* The relative power values are shown as a function of outlet 
or throat stagnation temperature for various values of stagnation pressure . 

These relative values are a combination of the relative energy ratios and the 
relative mass-flow ratios. The power values for the real gas are in general 
lower than those for the ideal gas. At the maximum pressure (fig. 9(d)) and 
minimum temperature, this reduction in the power required is about 9*5 percent 
for Mts = 0.2 (r = 1.025) and about 7.5 percent for Mxs =1-2 (r = 1.2). 

The shape of this power-ratio— temperature curve is essentially independent of 
the pressure ratio and/or Mach number. This independence is to be expected 
since the energy-ratio — temperature curve and the mass flow-temperature curve 
were essentially independent of the pressure ratio and Mach number, repectively. 
These power reductions due to real-gas effects are, of course, in addition to 
the large power reductions which result from operating at cryogenic temperature 
(fig. 1). 


APPROXIMATE METHODS 

For the engineering design of systems which utilize nitrogen, use of the 
complete calculation procedures of this report to include the real-gas effects 
on power calculations would be very cumbersome. Therefore, some approximate 
methods are now considered. 


Energy Per Unit Mass 

The following two equations are approximations for the energy required to 
compress a unit mass of real gas isentropically: 
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r " 


y 


r - 1 




(4) 


E 


r 


Z t,2 E i = Z t,2 




(5) 


Equation (4) is the energy per unit mass portion of the power equation given 
in reference 9. It is a result of the energy equation 

E r = J v dp 

and the assumptions that the pressure- temperature relationship for the ideal 
gas remains valid 

Y 

p = (Constant) 

and that Z varies linearly with pressure along the isentrope 
Z = (a + bp) g 


Equation (5) is derived with the same considerations, but Z is assumed 
to be constant along an isentrope. The data presented in figure 10 indicate 
that this is a reasonably good assumption. 

Reference 9 indicates the actual values of specific heat ratio y for the 
conditions prior to compression should be utilized in these equations. However, 
in reference 5, the isentropic expansion coefficients for nitrogen were found 
to remain near to the ideal diatomic gas value of 1.4. In addition, for this 
study, use of the ideal value of 1.4 gave more accurate results over the range 
of conditions considered herein than did use of the actual value of y . For 
equation (5), the use of either Z ^ ^ or Z^ g gives about the same degree 
of accuracy. For this report, Z^.^ -*- s use d’ 

Figure 11 shows a comparison of these approximate equations with the exact 
real-gas solutions. Both of the approximate solutions are within 0.5 percent 
of the exact values. Although equation (4) gives excellent values for the 
energy per unit mass, it is unnecessarily complex for the range of conditions 
considered herein. Simply multiplying the ideal values by the compressibility 
factor %t,2 gives results which are equally accurate. 
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Mass Flow 


The mass flow per unit area is m = pV. This equation may be rewritten 
in the following form: 


• ^ 

m = — PfMc 
p t 


The ratio of the real-gas mass flow 
given Mach number would be 


m„ 


to the ideal-gas mass flow m^ at a 


_ ( p/p t)r i_ 

m i ( p / p t)i Z t c i 


If p = p“ (Constant) is approximately valid for the real gas, and a remains 
near 1.4 (ref. 5), then 


” £ -ft 


Also 


( p / p t)r 

( p / p t)i 


1.0 (from ref. 5). With these considerations, 

•i 



The accuracy of this simple approximation is illustrated in figure 12. For the 
range of conditions considered in this report, this approximation is accurate 
to about 0.5 percent. 


Power 


r m r 1 r r , 

If — « Z f 9 and t— ^ ./ , then — «*' JZ*. 5 . Figure 13 illustrates 

E i t>2 m i p tf2 p i V t>2 

the accuracy of this simple approximation. For the range of conditions con- 
sidered herein,, the accuracy is within about 0.5 percent. 
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The results just presented indicate that simple corrections to the ideal- 
gas values of mass flow, compression energy, and compression power give accu- 
rate values for the real-gas case. 


CONCLUSIONS 

In this report, an analysis has been made of the real-gas effects of nitro- 
gen on the power required for steady-flow isentropic compressions. The compres- 
sions are assumed to occur at the fan of a transonic cryogenic wind tunnel. 

The analytic model tunnel was assumed to operate at stagnation temperatures 
from 310 K to saturation and at stagnation pressures from 1.0 to 8.8 atm. The 
results of this anlysis led to the following conclusions: 

1 . The energy to compress a unit mass of nitrogen at cryogenic tempera- 
tures is less than that required for the ideal gas. For the maximum pressure- 
minimum temperature conditions, this reduction in energy is on the order of 

14 to 17 percent. 

2. Tunnel mass flow at a given Mach number is higher for cryogenic nitro- 
gen than for the ideal gas. 

3. The power for isentropic compressions of cryogenic nitrogen is also 
less than that required for the ideal gas. At the maximum-pressure — minimum- 
temperature conditons, this reduction in power is on the order of 7.5 to 

9.5 percent. The power decrease is less than the energy decrease because of 
the increase in mass flow. 

A 

4. Simple compressibility factor corrections to the ideal values of mass 
flow, energy, and power give very accurate estimates of the real-gas values. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
January 18, 1977 
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APPENDIX A 


TUNNEL THROAT CONDITIONS - REAL-GAS CALCULATIONS 

Figure 14 contains a flow chart of the part of the program that makes 
the real-gas calculations of the throat conditions. The variables required 
for this part of the program are shown at the top of the chart. A step-by- 
step description of the calculations that occur at each block of the flow 
chart follows. All variables should be assumed to be real-gas variables 
except those specifically identified as ideal (subscript i) . 

Step 1 Subroutine PROP from the National Bureau of Standards (NBS) 

program is used to calculate the stagnation enthalpy hj. g and 

entropy S fc 2 , 

Step 2 An initial guess for the static conditions (p^, P^, and T^) is 

made by using the appropriate ideal-gas equations. 

Step 3 The real-gas value of static temperature T is initialized to the 
ideal value T^ and T is established as the iterative variable 
in order to force the solutions to converge on the desired M^. 

Step 4 Function DSFND iteratively solves for P from the entropy equa- 
tion Sf. ,2 = S( P,T). 

Step 5 Subroutine PROP gives values for static pressure p and static 

enthalpy h by using T and P from steps 3 and 4, respectively. 

Step 6 If the static temperature T and static pressure p are coincident 

with a point on the vapor pressure curve or if they lie in the liquid 
region, the solution is terminated. 

Step 7 Subroutine VSND (NBS program) calculates the velocity of sound c 
with T and P as inputs. 

Step 8 Subroutine MVCAL calculates the velocity and Mach number from these 
equations : 

V = \|2(h t , 2 " h) 



c 


Step 9 The mass-flow rate per unit area is now calculated 

m = PV 

Step 10 The calculated Mach number M is checked to see whether it is within 
0.00001 of the desired Mth* it is » the solution for the throat 

conditions is complete and the mass- flow rate is printed out. 
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APPENDIX A 


Step .11 If the Mach convergence is not satisfactory, subroutine TCHANG finds 
the slope of the temperature Mach number curve AT/AM for a constant 
entropy S fc 2 - A linear adjustment is made to the static tempera- 
ture T as’ required for Mach convergence. 

Step 12 This adjusted T is returned to step 3 for the next iteration, and 
steps 3 to 12 are repeated until the Mach convergence criterion is 
met at step 10. 


14 


APPENDIX B 


FAN CALCULATIONS - REAL GAS 

Presented in figure 15 is a flow chart of the part of the program that 
calculates the fan conditions for the real-gas (nitrogen) case. The objec- 
tives of this part of the program are to obtain an isentropic solution for 
the upstream stagnation quantities p^ i and Tj. ^ and then to calculate 
the energy per unit mass E and the pAwer P required for the compression. 

The downstream conditions p t 2 , ^ 2’ ^t 2’ an< * ^t 2* ** an P ressure 

ratio r, and the tunnel mass-flow rite m ’are required as inputs to this 
part of the program. The following is a step-by-step description of the cal- 
culations that occur at each block of the flow chart. Again, all variables 
should be assumed to be real-gas variables except those specifically identified 
as ideal (subscript i) . 

Step 1 The fan inlet pressure is calculated directly as shown. Note that 

Pj. ■) and p t 2 have the same values for both the real- and ideal- 

gas cases. ’ 

Step 2 An initial guess for the fan inlet temperature (T^ is made 

using the ideal-gas equation. ’ 

Step 3 T^ i is set up for iteration in order to force the upstream 

entropy S t , to be the same as the downstream entropy S t 2 (see 

step 6). ’ 

Step 4 Subroutine PROP is called with p t 1 and T t ■, to give values for 
the upstream enthalpy h^. ^ and entropy S t j . 

Step 5 The upstream temperature is increased by a factor of 1.001 and is 
used in subroutine PROP with p^. ^ to give another point of the 
entropy- temperature curve ^tjl'j'u' 

Step 6 A comparison of S t ^ and S t 2 is made to see whether these 
two entropies agree’ to within 1 * 10-8. 

Step 7 If the two entropies have not converged sufficiently, a temperature- 
entropy slope is calculated from the information of steps 4 and 5. 

Step 8 This temperature- entropy slope is used to adjust T fc ^ to account 
for the difference in the inlet and outlet entropies’ S^ >2 " St ? -|, 
and the adjusted T^i is returned to step 3 for the next iteration. 

Step 9 When the convergence criteria of step 6 are finally satisfied, the 
energy and power are calculated. 

Step 10 The solution is complete, and the fan inlet conditions and power 
values are stored for utilization in other parts of the program. 
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APPENDIX C 


PROGRAM LISTING AND SAMPLE OUTPUT 

Presented here is a listing of the program that was developed in order 
to calculate the power and energy required to compress nitrogen and the ideal 
diatomic gas isentropically. This program is a modification of those of refer- 
ences 5 and 7. In addition, a sample output is presented. 


Program Listing 


PROGRAM N2P0WER ( INPUT » OUTPUT *TAPE5=INPUT»TAPE3=0UT PUT) 

CALL DATA N2 

THIS PROGRAM CALCULATES THE POWER REQUIREMENTS FOR ISENTROPIC 
COMPRESSIONS OF BOTH NITROGEN AND AN IDEAL DIATOMIC GAS 

REAO 100, NC 
100 FORMAT (15) 

C THESE POWER REQUIREMENTS ARE CALCULATED FOR A SPECIFIED NUMBER OF 
C CASES (NC)* EACH OF WHICH INVOLVES KNOWN VALUES FOR THE TEST-SECTION 
C MACH NUMBER (BM2) * FAN OUTLET PRESSURE (POT)* AND THE RATIO OF FAN 
C OUTLET TO INLET PRESSURE (PR4) 

C 

C THE POWER AND MASS-FLOW RATE CALCULATIONS INVOLVE CROSS-SECT I ONAL 
C AREAS WHICH CORRESPOND TO THE THROAT AREA OF THE TEST SECTION f 

C 

DO 500 L = 1 * NC 
REAO 200* POT* PR4* BM2 
200 FORMAT (3F10.0) 

PRINT 210* BMP, -PR4 * POT* (P0T/PR4) 

210 FORMAT < 1 HI . 32X **ENEPGY AND POWER REQUIREMENTS FOR FAN-DRIVEN CRYO 
1GENIC WIND TUNNEL */49X**ISENTR0PIC COMPRESSION IN NITR0GEN*///2X** 
2 YE ST- SECT I ON MACH NUMPER =* * F5 . 2/?X * *F AN PRESSURE RAT IO** 7X »*=* * F7 
3.4** */2X**FAN PRFSSURE. 0UTLET**5X ,*=»,E7.4,* ATM*/2X**FAN PRE 

4SSURE* INLET*,6X»*=**F7.4,* ATM*///2X»*E = FAN ENERGY*/2X » *P = FAN 
5 P0WER*///11X*17(*-*) **REAL GAS* . 1 7 (*-*)» 4X . 1 7 (*-*). *IDEAL GAS*»16 
6 (*-*) »4X .5 (*-*) **REAL/ IDEAL VALUES* *4 ( *-*) /2X « *TT * IN**4X * 2 (*TT »OUT 
7*.6X**MD0T**7X»*E/MASS**6X,*P/AREA*.5X) * 1X»*MD0T*»5X**E/MASS*»4X»* 
flO/AREA*/4X.*K**8X*2 (*K**7X**KGM/S-M2**5X»*J/KGM**7X.*J/S-M2*»7X) /)_ 

C 

C FOR EACH CASE* THE CALCULATIONS ARE PERFORMED FOR EACH OF 
C INCREMENTALLY DECREASING VALUES OF FAN OUTLET TEMPERATURE (TOT) UNTIL 
C SATURATION CONDITIONS ARE REACHED 
C 

DO 280 1=1.50 
IF ( I .GT. 20) GO TO 205 
TOT= 320. - I * 10 
GO TO 206 

205 TOT = 120. - (1-20) * 5 


16 



onooooooooo o o o oooooonoo onoo oooooooooo 


APPENDIX C 


THE LIQUID REGION CONDITION IS REACHED WHEN THE TEMPERATURE FALLS 
BELOW THAT OF THE CRITICAL POINT* AND THE PRESSURE RISES ABOVE THE 
CORRESPONDING SATURATION PRESSURE (PSAT) 

THE CRITICAL POINT FOP NITROGEN IS AT PRESSURE = 33.5 ATM AND AT 
TEMPERATURE = 126 DEGREES KELVIN 

WHEN THIS SITUATION OCCURS* THEN THE NEXT CASE MUST BE CONSIDERED 
206 PSAT =33,5 

IF (TOT .LT. 126) PSAT = VPN(TOT) 

IF (POT .GT. PSAT) GO TO 385 

USING THE KNOWN PROPERTIES OF NITROGEN* THE DENSITY (DOT). ENTHALPY 
(HOT)* AND ENTROPY (SOT)* AT THE OUTLET OF THE FAN. CAN BE CALCULATED 

CALL PROP (TOT»POT*DOT,l.HOT*SOT*U) 

STATIC CONDITIONS AT THE TUNNEL THROAT WHICH CORRESPOND TO THE DESIRED 
THROAT MACH NUMBER ARE ESTIMATED BY ASSUMING ISENTROPIC FLOW OF IDEAL 
GAS FROM FAN OUTLET TO THROAT 

T IS THE TEMPERATURE 
D IS THE DENSITY 
P IS THE PRESSURE 

BM3=BM2 

IF (BM3.GT.1.0) BM2= 1.0 
T n 1 , / ( 1 . *BM2**2. /5. ) *TOT 
D = DOT*(T/TOT> **2.5 
P = POT* (T /TOT) **3.5 

ANOTHER TEST FOR SATURATION CONDITIONS 

PSAT = 33.5 

IF (T.LT.126.) PSAT = VPN(T) 

IF (P .GT, PSAT) GO TO 385 

THE THROAT TEMPERATURE IS VARIED UNTIL THE REAL- GAS ISENTROPIC 
SOLUTION CONVERGES ON THE DESIRED THROAT M ACH NUMBER 

H IS THE ENTHALPY 

S IS THE ENTROPY 

BMACH IS THE TRIAL MACH NUMBER 

VEL IS THE FLOW VELOCITY 

QN2 IS THE MASS-FLOW RATE FOR THE REAL-GAS CASE 

DO 217 K = 1*100 
D = DSFND(SOT*T.O) 

CALL PR0P(T,P,D.2» H.S.U) 
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ANOTHER TEST FOR SATURATION CONDITIONS 
0SAT = 33.5 

IF (T.LT. 126. ) PSAT = VPN(T) 

IF (P .GT. PSAT) GO TO 385 
CALL VSNDCT.P.D.2.W) 

CALL MQRCAL <H0T.H,D.T,W,BMACH,QN2.REY) 

call mvcalihot.h.w.bmach.vel) 

QN2= 28.0134 *D*VEL 

IF CABS(BMACH-BM2).-. 00001) 218, 218. 216 

216 IF (K.EQ.100) GO TO 330 
CALL TCHANG (BMACH.BM2.T.D.S0T.H0T) 

217 CONTINUE 

THESE ARE THE FAN INLET PRESSURE (P) AND AN ESTIMATE OF THE FAN 
INLET TEMPERATURE (T) 

218 P = P0T/PR4 
T = T0T/PR4**.2857 

NOW THE PROGRAM VARIES FAN INLET TEMPERATURE SO THAT THE INLET 
ENTROPY CONVERGES TO THE OUTLET AND/OR TEST-SECTION VALUE (ISENTROPIC 
COMPRESSION) 

TIJ IS THE INCREMENTED TEMPERATURE 
SU IS THE CORRESPONDINGLY INCREMENTED ENTROPY 

DO 260 J=1.100 
TU = 1.001 * T 
CALL PROP (T.P.D.l.H.S.U) 

CALL PROP (TU.P.DU. 1 .HU.SU.U) 

IF (ABS(S-SOT) -l.E-08 * SOT) 265.265.256 
256 T= T ♦ (TU-T)/(SU-S) * (SOT-S) • .6 
IF ( J .EQ. 100) GO TO 406 
260 CONTINUE 

NOW THAT ALL OF THE NECESSARY CONDITIONS FOR THE TUNNEL HAVE BEEN 
DETERMINED. THE POWER REQUIREMENTS CAN BE CALCULATED AND PRINTED OUT 

EN1 IS THE ENERGY REQUIREMENT PER UNIT MASS FOR THE REAL-GAS CASE 
265 EN1 = 35.697* ( HOT-H ) 

C R IS THE NITROGEN GAS CONSTANT. IN J/KGM-K 
R = 296.791 

C TOI IS THE OUTLET TEMPERATURE EOR THE IDEAL-GAS CASE 
TO I = PR4**.2857*T 

C EN2 IS THE ENERGY REQUIREMENT PER UNIT MASS FOR THE IDEAL-GAS CASE 
EN2 = 3. 5*R*T* (PR4**. 2857-1 .0) 

C EN3 IS THE RATIO OF REAL TO IDEAL ENERGY PER UNIT MASS REQUIREMENT 
EN3 = EN1/EN2 
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C HP 1 IS THE FAN POWER REQUIRED PER UNIT CROSS-SECTIONAL AREA FOR THE 

C REAL-GAS CASE 

HP 1 = EN1*QN2 

C QN2I IS THE MASS-FLOW RATE FOR THE IDEAL CASE 

PN2I = SQRT(1.4/<9.74015E-11*R*TOI) ) *P0T*BM2/ ( 1 . ♦BM2**2 . /5 . 0 ) **3. 
PM2=BM3 

C QN3 IS THE RATIO OF REAL TO IDEAL MASS-FLOW RATE 
QN3 = QN2/QN2 I 

C HP2 IS THE FAN POWER REQUIRED PER UNIT CROSS-SECTIONAL AREA FOR THE 

C IDEAL-GAS CASE 

HP2 = EN2*QN2I 

C HP3 IS THE RATIO OF REAL TO IDEAL POWER REQUIRED PER UNIT CROSS- 
C SECTIONAL AREA 

HP3 = HP1/HP2 
C 

PRINT 2TO«TfTOT*QN2*ENltHPl*TOI»QN2I*EN2»HP2tQN3»EN3»HP3 
27 0 FOPMAT (F7.1 » 2 ( F9 . 1 * 2F 1 2 « 1 * E 1 3 , 4 ) ,3F10.4) 

280 CONTINUE 

SO TO 500 
330 °RINT 340 

340 FORMAT ( * TSMP-MACH DID NOT CONVERGE IN LOOP 217*) 

GO TO 500 
385 PRINT 390 

390 FORMAT<* SATURATION CONDITIONS REACHED AT THROAT MACH NUMBER*) 

GO TO 500 

406 PRINT 407 

407 F ORM AT ( * TEMP-ENTROPY DID NOT CONVERGE IN LOOP 260*) 

500 CONTINUE 

STOP 

END 


SUBROUTINE CPVTD ( T » D . CP » C V ) 

COMMON /CRPR/ CR ( 3 ) /METH/ M 
COMMON /RFPR/ RF(10) 

C 

C.... ROUTINE TO CALCULATE C V AND CP FROM THE EQUATION OF STATE 
C..., WRITTEN 7/21/71 A MYERS 
C.... REVISED 12/5/71 A MYERS 
C 

DC = CR <2) 

TC=CR (3) 

R =RF ( 5 ) 

AK=RF (6) 

CALL PFND(T»D»P) 

TF(M.EQ.l)GO TO 1 
IF (T.GT.TC) GO TO 1 
IF (D.GT.DC)GO TO 2 
1 CVO=CPIG(T)-(R*AK) 

CV=CVO- (FING3 (T»D) “FING3 (T»0.00) ) * AK 
3 F 1 =T/D**2 
F2=DPDT(T.D) 

F 3 = DPDD ( T * D ) 

CP = CV» <F1*F2**2/F3)*AK 
RETURN 
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ht=o.i 

T1=T*DT 


T2=T-DT 

CALL LPR0P(T1,P1,D.1.H,S.U1> 
CALL LPR0P(T2»P2»D»1»H»S»U2) 
CV=(U1-U2)/(2.00*DT) 

GO TO 3 
END 


SUBROUTINE DAT AN2 

COMMON /CEOS/G (41) /CVPN/GV(11) /CIGCP/GI (9) /CSL/CL<7) /CSV/CV(7) 

. /CRPR/CR ( 3 ) /CTEVP/CT ( 8 ) /RFPR/RFUO) /METH/ M 

.... IF THE PROPERTIES OF NITROGEN ARE TO BE calculated^, A CALL TO this 
.... SUBROUTINE MUST BE THE FIRST CALL STATEMENT IN THE MAIN PR06RAM 

.... THE COMMON BLOCKS INITIALIZED IN THIS ROUTINE HOLD THE FOLLOWING 
.... INFORMATION - 

/CEOS/ G ( 4 1 ) COEFFICIENTS OF THE EQUATION OF STATE 

/CIGCP/ Cl G ( 9 ) COEFFICIENTS OF THE IDEAL-GAS HEAT CAPACITY EQUATION 


/C VPN/ GG(ll) COEFFICIENTS OF THE VAPOR PRESSURE EQUATION 


/CRPR/ CR ( 3 ) THE CRITICAL PROPERTIES IN THE SAME UNITS AS THE 
EQUATION OF STATE * 

CR ( 1 ) “CRITICAL PRESSURE 
CR(2)=CRITICAL DENSITY 
CR (3) “CRITICAL TEMPERATURE 

/CSL/ SLC7) COEFFICIENTS OF EQUATION USED TO APPROXIMATE 
THE SATURATED LIQUID DENSITY AS A FUNCTION OF 
TEMPERATURE 


/RFPR/ RF(10) REFERENCE PROPERTIES 

RF ( 1 ) “REFERENCE ENTHALPY AT TEMPERATURE TOH 
RF(2)=RFFERENCE ENTROPY AT TEMPERATURE TOS 
RF (3) *TEMPERATURE TOH 
RF (4) ^TEMPERATURE TOS . 

RF(5)=GAS CONSTANT IN UNITS OF EQUATION OF STATE - R 
RF (6 ) sCONVERSION FACTOR TO CHANGE UNITS OF 

THE EQUATION OF STATE TO DESIRED ENERGY UNITS 
RF ( 7 ) =MOLECUL AR WEIGHT . 

RF(8)=TRIPLE POINT TEMPERATURE 
RF ( 9 ) - NOT USED 

RF (10) - NOT USED 

/METH/ M INDICATES METHOD TO BE USED IN THE CALCULATION 

OF LIQUID PROPERTIES 

M= 1 INDICATES ISOTHERM INTEGRATION THROUGH THE DOME 
M=2 INDICATES THE USE OF THE CLAPEYRON EQUATION 
THROUGH THE DOME 
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Msl 

G( 1)= 0.1362247692728270-02 

G ( 2>= 0*1070324699085910 00 

G ( 3 ) a -0.2439007218714130 01 

G( 4)* 0.341007449376470D 02 

G( 5)= -0.4223743094661670 04 

G ( 6)= 0.1050986002464940-03 

G ( 7) = -0,1125948265220810-01 

G ( 8 ) = 0.1426007892709070-03 

G< 9)= 0.1846905016090070 05 

G ( 1 0 ) = 0.8111400825887760-07 

0 ( 1 1 ) e 0.2330116450380060-02 

G ( 12) = -0.5077525863509860 00 

G ( 1 3 ) a 0.4850278819312140-04 

G ( 14 ) s -0,1136567641153640-02 

G ( 15) = -0.7074302735405750 00 

G ( 1 6 ) 3 0.7517066488526800-04 

G ( 1 7 ) 3 -0.1116141195374240-05 

G ( 18 ) = 0.368796562233495D-03 

G ( 1 9 ) = -0. 20 f3 176913477290- 05 

G ( 20 ) = -0.16971 744 475594 9D 05 

G ( 21 ) = -0.1 1 97 19240 044 192D 06 

G ( 22) = -0.9752182720382010 02 

G ( 23 ) = 0.55463971315182 3D 05 

G ( 24 ) * -0.1799204504434700 00 

G ( 25 ) = -0.2565829260771840 01 

G ( 26 ) = -0.4137077150907890-03 

G ( 27 ) s -0.2562454153002930 00 

G ( 28 ) = -0.12422237 37 4006 3D-06 

G ( 29 ) 3 0.1035565358401650-04 

G ( 30 ) = -0.5386991665583030-09 

G (31 ) = -0,7574154 1 2839596D-08 

G ( 32 ) = 0.5853671720695210-07 

G ( 4 1 ) 3 -0.5600000000000000-02 

GV(1>3 0.83944094440 04 
GV(2)3-o. 18900452590 04 
GV (3)3-0.72822291650 01 
GV ( 4 ) =0 • 0D0 

GV ( 5 ) = 0.55560638250-03 
GV (6) s-o.59445446620-05 
G V ( 7 ) a 0.27154339320-07 
GV (8) 3-0.48795359040-10 
G V (9)3 0.5095360824D 03 
GV(10)3 0.10228509660-01 
GV(11)3 0.1950000000000000 01 

PI(1)3 -0.7352104011572520 03 

G I (2)3 0.3422399804119780 02 

G 1 (3)= -0.557648284567620D 00 

G I (4 ) * 0.3504042283087560 01 

G I (5)= -0.1733901850810050-04 

GI (6)3 0.1746508497664630-07 

GI (7)= -0.3568920335443480-11 

Gl (8)3 0.1005387228088340 01 

G I (9)s 0 . 33534061 00 000 OOD 04 

CL ( 1 ) 3 0.1942440319220000 02 

CL ( 2 ) 3 0.5708374819420000 02 

CL ( 3 ) = -0.2432698504630000 03 


L 
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CL (4) = 

0.885168381502000D 

03 

CL (5) = 

«0. 1639367977 34 000D 

04 

CL (6) = 

0,11574 32 00533000D 

04 

CL (7) * * 

0.1 0 182209832700 OD 

01 

CV(1)= 

0. 1 633334553880 00D 

01 

CV(2)= 

-0.940437709170000D 

01 

CV(3>= 

0. 21 8527460 266 000 D 

02 

CV<4)= 

-0.102687437637000D 

03 

CV(5)= 

0. 1 879497448620 00D 

03 

CV (6) = 

-0 • 164024367971 00 OD 

03 

CV (7) = 

-0. 9573.1639031 6000D-01 

CR ( 1 ) =33. 

555D0 


CR(2)=11. 

2100 


CR (3 ) = 1 26 

• 200 


CT(1>= 

-0. 14206478620000 0D-02 

CT(2)= 

0.1 29088086900 000D-01 

CT ( 3 ) = 

0.0 


CT (4) = 

0.0 


CT (5) * 

0.0 


CT (6) = 

0.0 


CT ( 7 ) = 

0.0 


CT (8) = 

0.0 


RF ( 1 ) =8.669D03 


RE(2)=191 

.502000 


RF (3) =298 

.15000 


RF(4)=298 

.15000 


RF (5 ) * 

0.820539000000000D- 

■01 

RF (6) = 

0.1013278000000000 

03 

RF<7)=28. 

0134D0 


RF ( 8 ) =63. 

1 4800 


RETURN 
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SUBROUTINE DCALC (D » • . P » DL ♦ DH ) 

DATA (MAX=30) 

DATA(EPS=l.E-5) 

C 

C.... ROUTINE TO PERFORM ITERATIVE SOLUTION OF THE EQUATION OF STATE 
C 

C.... DL IS LOWER BOUND ON DENSITY 
C.... DH IS UPPER ROUND ON DENSITY 
C 

C.... THE DESIRED DENSITY MUST LIE BETWEEN DL AND DH 
C 

C.... ALGORITHM IS MODIFIED VERSION OF 
C 

C.... * A QUADRATIC FORMULA FOR FINDING THE ROOT OF AN EQUATION* BY 

C.... LI. G. CHAMBERS MATHEMATICS OF COMPUTATION VOL 25 NO 114 APRIL (1971) 

C 

C.... WRITTEN 2/7/72 A MYERS 
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IDS = 0 
IC = 0 

2 CONTINUE 
0 1 =DL 
D2=DH 

CALL PEND (T.D1 »Pl ) 

CALL PEND ( T . D2 • P2 ) 
y l =p I -p 

Y2=P2-P 
1 IC=IC*1 

IF (IC.GT.MAX)GO TO 5 
t F ( I DS.EQ » 1 ) GO TO 7 
IF ( IC . GT • ?0 ) GO TO 6 
ns= (Dl*Y2-D2*Yl > / (Y2-Y1 > 

GO TO 8 

7 DS= (D1*D2> /2.00 
fi CONTINUE 

CALL PFND(T.DS.PS) 

YS=PS-p 

03=DS*Y1*Y2/ ( (YS-Y2) * (YS-Y1 > ) ♦ 

. D1*YS*Y?/ ( (Y1-Y2) * ( Yl-YS) > ♦ 

. D2*YS*Yl/( (Y2-Y1 > * ( Y2-YS) > 

CALL PFND (T.D3*P3) 

Y3=P3-P 

I F ( ABS(Y3) ,LE.EPS)GO TO 3 
IF (Y3.GT. 0.00)60 TO 12 
IF (Y3.LT.YDGO TO 12 
Y 1 = Y3 

01 = D3 

12 TF(YS.GT. 0.00)60 TO 13 
IF (YS.LT »Yl ) 60 TO 13 

Y 1 = YS 
ODDS 

13 TF (Y3.LT. 0.00)60 TO 14 
IF ( Y3.GT .Y2) GO TO 14 
Y2 = Y3 

D2 = D3 

14 IF (YS.LT. 0.00)60 TO 1 
IF (YS.GT.Y2) GO TO 1 
Y? = YS 

02 = DS 
GO TO 1 

3 CONTINUE 
0 = D3 
RETURN 

5 WRITE (3»300) T»P»DL«DH.D1 »Y1,D2.Y2»DS»YS»D3*Y3 
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300 F ORMAT ( * DCALC FAILED TO CONVERGE 



« 

T 

= 

**F15.7*/t 


* 

P 

= 

*»F15.7*/ t 



DL 

= 

*»F15.7»/« 


. * 

DH 

= 

*»F15.7»/, 


« 

D1 

= 

*»F15.7»/ . 



Yl 

= 

*,F15.7./, 


« 

D2 

= 

*»F15.7»/» 


* 

Y? 

= 

*»F15.7»/» 


» 

DS 

= 

*»F15.7*/» 


* 

YS 

s 

*»F15.7*/» 


* 

D3 

= 

*»F15.7»/* 



Y 3 

= 

*,F15.7) 

D ETU»N 





CONTINUE 




ms=i 

IC = 0 
GO TO 
END 

2 





*/ 


SUBROUTINE DFND ( T » P * D ♦ K ) 

COMMON /RFPR/ RF (10) 

COMMON /CRPR/CP ( 3 ) 

.... ROUTINE TO GENERATE TRIAL DENSITIES FOR ITERATIVE SOLUTION OF 
.... THF EQUATION OF STATE FOR DENSITY GIVEN TEMPERATURE AND PRESSURE 

K =0 INDICATES SINGLE PHASE POINT 

K =1 INDICATES T ♦ P ARE FOR THE SATURATED LIQUID 
K r? INDICATES T ♦ P ARE FOR THE SATURATED VAPOR 

.... WRITTEN 2/10/72 A MYFRS 

IF ( (K.LE.2) .AND. (K.GE.O) )GO TO 1 
WR I TE ( 3 » 30 0 ) K 

300 F ORMAT ( 27H *** ERPOR IN CALL DFND ****/f 
. * K MUST EQUAL 0» 1* OR 2*»/t 

. * K = *»I10) 

RFTURN 
1 PC=CR(1) 

DC=CR (2) 

T C=CR ( 3 ) 

IF(K.GT.O)GO TO 5 
IF (T.GE.TC)GO TO 2 
VP=VPN(T) 

T F ( P • LE • VP ) GO TO 3 
4 DH=3.100*DC 
0L=DSATL(T> 

CALL DCALC (D,T.P»DL.DH) 

RFTURN 
3 OL = 0 

DH=DSATV (T) 

CALL DCALC (D.T»P*DLtDH) 

RETURN 
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2 OL = 0 

DH=3.100*DC 

IF ( ( T .6T .1000.00) .AND. (P.LT.300.00) >DH=DC 
CALL DCALC(D.T»P.DL.DH> 

RETURN 

5 IF (K .EO. 1 ) 60 TO 4 
GO TO 3 
END 


SUBROUTINE DPDTVP (T.P.DPDT) 

COMMON/CVPN/G<ll> /CRPR/CP (3) 

COMMON /SCRH/ X(40) 

.... CALCULATE DP/DT FOP THE VAPOR PRESSURE EQUATION 

.... WRITTEN 7/22/71 A MYERS 

TC=CR<3) 

A=G(11) 

T?»T*T 
T3=T*T2 
T4=T*T3 
TS=T*T4 
X ( 1 ) =-l , 00/T2 
X ( 2) =0 
X < 3 ) =1 . 00 
X (4) =2.00*T 
X(5)=3.00*T2 
X <f>) =4.00*T3 
X (7)=5.00*T4 
X (8) =6.00*T5 
X (9) =1 . 00/T 

X (10)s(TC-T)**(A-1.00)*(-A) 

OPDT = 0 
DO 1 1=1*10 
1 OPDT=DPDT ♦X(I)*G{I) 

0R0T=DPDT*P 

RETURN 

END 


SUBROUTINE LPROP (T.P.D.KtH.S.U) 

COMMON /RFPR/ RF(10) 

COMMON /METH/ M 

. ... ROUTINE TO CALCULATE THE PROPERTIES OF THE LIQUID 
. ... WRITTEN 7/20/71 A MYERS 
.... REVISED 12/8/71 A MYERS 

K =1 INPUT IS T ♦ 0 
K =2 INPUT IS T ♦ P 
K =3 INPUT IS T. P* ♦ D 

TF(M.EQ.2)G0 TO 1 

CALL VPROP<T,P,D.kVH,S»U> 

RETURN 
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1 R =PF (5) 

AK=RF(6) 

IF(K.E0.1)CALL PFND ( T ♦ D t P ) 

I F ( K . EQ • 2 ) CALL DFND ( T • P • D » 0 ) 

V D = VPN ( T ) 

CALL DFND(T»VP»DSV»2) 

CALL VPR0P(T»VPP»DSV*1»HSV»SSV*USV) 

VSV=1 • OO/DSV 

CALL DFND(T,VP*DSL.1> 

VSL=1. 00/DSL 
CALL DPDTVP(T,VP.DPDT) 
HSL=HSV-T*DPDT* (VSV-VSL)*AK 
SSL=SSV* (HSL-HSV) /T 

USL=USV* (HSL-HSV) - (VP* (VSL-VSV) ) *AK 
OLD=ALOG(D) 

DLS = ALOG(DSL) 

F1D=R*DLD-FING1 ( T ♦ D ) 

F1S=R*DLS-FING1 ( T » DSL ) 

F?D=FIHG2(T,D) *R*T*DLD 

F2S=FING2 (T»DSL) ♦R*T*DLS 

S=SSL-(F1D-F1S)*AK 

U = USL* ( (F2D-F2S)-T*(F1D-F1S) )*AK 

H=U ♦ (P/D ) * AK 

RETURN 

FMD 


SUBROUTINE MQRCAL ( HOT * H » D » T »W»BMACH»GN2»REY) 

SUBROUTINE CAL M.Q.REY/M 

VEL = SORT (2/. 0280134* (HOT-H) ) 

BMACH = VEL/W 

QN2 = D*28.0134*VEL**2/2.0 

DD=D*. 0280134 

EMU =VISC (DD*T) *1 .E-04 

REY = 28.0134*D*Bm'ACH*W/EMU 

RETURN 

END 


SUBROUTINE MVCAL ( HOT » H * W * PMACH * VEL ) 
VEL = SORT (2/. 0280134* (HOT-H) ) 

BMACH = VEL/W 

RETURN 

END 


SUBROUTINE PFND(T»D»P> 
COMMON /RFPR/ RF(10) 
COMMON /CEOS/G (41) 
COMMON /SCRH/ B ( 40 ) 


non 
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.... THIS ROUTINE CALCULATES PPESSURE GIVEN TEMPERATURE AND DENSITY 
.... FROM THE EQUATION OF STATE 

R=RF (5) 

02 = 0*0 
D3=D2*D 
04=03*0 
05=D4*D 
06=D5*D 
07=D6*D 
D«=D7*0 
D9=D8*D 
D10=D9*D 
01 1=D10*D 
012 = 011*0 
Q13=D12*D 
TS=SQRT (T) 

T2=T*T 
T3=T2*T 
T4=T3*T 
GM=G(41 ) 

F = EXP (CM *D2 ) 

8( 1 ) =D2*T 
9( 2) =D2*TS 
B( 3) =D? 

9( 4 ) =D2/T 
9( 5>=D2/T2 
R( 6 ) =D3*T 
9( 71=03 
9( 8>=D3/T 
9( 9 1 =D3/T2 
9(10) =D4*T 
9(11) =D4 
9(12) =D4/T 
9(13) =05 f 
9(14) =D6/T ' 

9(15) =D6/T2 
9(16) =D7/T 
9(17) =D8/T 
9(18) =Q8/T2 
9(19) =D9/T2 
9(20) =D3*F/T2 
9(21) =D3*F/T3 
9(22) =05 *F /T? 

9(23) =05*F/T4 
9 ( 24 ) =07*F/T2 
8 (25) =D7*F/T3 
9(26) =09*F /T2 
8(27) =09*F/T4 
B(28)=DU*F/T? 

8 (29) =01 1*F/T3 
9(30) =D13*F/T2 
8(31) =01 3*F/T3 
9(32) =01 3*F/T4 
N=32 
P=0 

DO 1 1 = 1 .N 
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1 

S 


1 p=P*B(I)*G(I) 
p=P+r*D*T 
RETURN 
END 


SUBROUTINE PROP (T.P.D.K.H.S.U) 

COMMON /CRPR/ CR ( 3 ) /METH/M 
C 

C.... GENERALIZED PROPERTY CALCULATOR 
C.,.. WRITTEN 7/20/71 A MYERS 
C.... REVISED 12/10/71 A MYERS 
C 

C.... ROUTINE CALCULATES PROPERTIES FOR FOLLOWING INPUT OF K 
C 

C K = 1 INPUT IS T ♦ P RFTURNS D, H, S, ♦ U 

C K=2 INPUT IS T * D RETURNS P, H, S. ♦ U 

C K = 3 INPUT IS T RETURNS P, D» H, S. • U FOR SATURATED VAPOR 

C K = 4 INPUT IS T RETURNS P. D. H ♦ S * • U FOR SATURATED LIQUID 

C 

C.... NOTE: ALL REAL VARIABLES IN CALL STATEMENTS TO ROUTINES IN 

C THIS PACKAGE MUST BE TYPE REALMS (DOUBLE PRECISION) 

C 

C.... NOTE: THE FIRST CALL STATEMENT IN THE USER*S MAIN PROGRAM MUST 

C BE TO A DATA INITIALIZATION ROUTINE 

C EXAMPLE: CALL DAT AN2 

C 

C.... NOTE; THE METHOD OF PROPERTY CALCULATION IS DETERMINED BY THE 
C VALUE OF M CONTAINED IN COMMON BLOCK /METH/M 

C 

C M=1 INDICATES PROPERTY CALCULATION TO BE CARRIED OUT BY 

C CONTINUOUS INTEGRATION OF ISOTHERMS THROUGH THE TWO PHASE 

C REGION 

C M=£ INDICATES PROPERTY CALCULATION IS INTERRUPTED AT THE TWO 

C PHASE VAPOR BOUNDARY AND THE CLAPEYRON RELATION WITH THE 

C VAPOR PRESSURE EQUATION IS USED TO CALCULATE THE LATENT 

C HEAT. INTEGRATION OF ISOTHERMS IS CONTINUED AT THE 

C SATURATED LIQUID BOUNDARY 

C 

PC=CP(1) 

OC=CR ( 2 ) 

TC=CR(3> 

IF ( (K.GT.O) .AND. (K.LT.5) ) GO TO 1 
WRITE ( 3 * 30 0 ) K 

300 F ORMAT ( 27H EPPOR IN CALL PPOP ***,/» 

. * K MUST FGUAL 1.2.3. OR 4*,/, 

. * K = *.110) 

RETURN 

1 IF ( K .LT • 3 ) GO TO 3 
IF(T.LE.TC)GO TO 2 
WRITE (3,301 ) T 

301 FORMAT (27H' *** FPROR IN CALL PROP ***»/» 

. * SATURATION PROPERTIES HAVE BEEN REQUESTED*./, 

. * FOP A TEMPERATURE THAT EXCEEDS CRITICAL*./, 

. * t = *»F15.5) 
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2 P = VPN ( T ) 

IF(K.E0.3)CALL DFMD (T»P.P.2) 
IF (K.EQ.4)CALL DFNO (TtP.D. 1) 
GO TO 4 

3 IF (K«GT« 1 ) GO TO 7 
CALL QFND<T«P»0»0) 

IF (T.GT.TC) GO TO 5 

4 IF ( D • GT « DC ) GO TO 6 

5 CALL VPROP (T»P»D«3»H»StU) 

return 

8 TF(M.E0.1)G0 TO 5 

CALL L p R0P(T»P«0«3*H»SfU) 
RETURN 

7 IF(T.GT.TC)GO TO 8 
VPsVPN(T) 

CALL DFNO <T tVP»DV»2) 

CALL DFND(T,VP*OL»l) 

IF (D.GF.DL) GO TO 9 
IF (D.GT.DV) GO TO 10 

8 CALL VPROP (T*P«0. 1 «H«SfU) 
RETURN 

9 IF (M«EQ» 1 ) GO TO 8 

CALL LPROR (TtPtO, 1 »H»S.U) 
RETURN 

10 VL=1 • OO/OL 
VV=1 .00/DV 
V =1.00/0 
X=(V-VL)/(VV-VL> 

CALL VPROP (T » P « DV » 1 »HV»SV«UV) 
IF (M.EQ.2) GO TO 11 
CALL VPROP(T.P»OL»l.HL.SL*UD 
GO TO 12 

11 CALL LPROP(T*P,OL*l .HL»SL»UL) 

12 M=HL*X* (HV-HL) 

S=SL*X* (SV-SL) 

U=UL*X* (UV-UL) 

RETURN 

END 


SURROUT INE TVP(P.T) 

COMMON /CTEVP/GT ( P ) /CRPP/CR(3) 

COMMON /SCRH/ X ( 4 0 ) 

ROUTINE TO SOLVE VAPOR PRESSURE EQUATION ITERATIVELY FOR 
TEMPERATURE BY NEWTON*S METHOD 
WRITTEN 7/22/71 A MYERS 


TC=CR(3) 
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.... USE TEMP EXPLICIT EON FOR FIRST APPROX 

P2=P*P 

P3=P2*P 

P4=P3*P 

P5=P4*P 

P6=P5*P 

X(l)=ALOG(P) 

X (2) =1 ,00 
X(3)=P 
X(4)=P2 
X (5) =P3 
X(6)=P4 
X (7)=P5 
X(8)=P6 
T = 0 

00 1 1 = 1.8 

1 T=T*X ( I ) *GT ( I ) 

T= 1 . 00/T 

.... T IS NOW FIRST EST OF T 

I TRMAX=25 
EPS=l.E-7 

DO 2 ITER=1 . ITRMAX 
PPsyPN ( T ) 

CALL DPDTVP(T.P.DPDT) 

OFLTA=(P-PP)/DPOT 
T = T ♦DELT A 

I F ( ABS(DELTA/T) .LT .EPS > RETURN 

2 CONTINUE 

WRITE (3»300)P.T. DELTA 

300 FORMAT (25H *** TVP DID NOT CONVERGE,/* 
. * P =*,F15.7, 

, * T =*,F15.7, 

. * DEL =* , F 1 5 . 7 ) 

RETURN 

END 


SUBROUTINE TCHANG ( RMACH , AM7 , T , D . SOT , HOT ) 

C THIS SUB, CHANGES TEMP ACCORDING TO TEMP VS MACH CURVE 

TUP =T*.00001*T 
TDN=T-.00001*T 
DijP=DSFND(SOT»TUP,D) 

D0N=DSFND(S0T,TDN,DUP> 

CALL PROP (TUP. PUP, DUP, 2, HUP, SUP .UUP) 

CALL VSND(TUP,PUP,DuP,2,WUP) 

CALL MVCAL (HOT ♦ HUP ♦ WUP » UMACH » UVEL ) 

CALL PROP(TDN,PDN,DDN,2,HDN,SDN,UDN) 

CALL VSND (TON,PDN,DDN,2,WDN) 

CALL MVCAL (HOT, HDN, WON, DMACH.DVEL) 

T=T- (8MACH-AM7) * (TUP-TDN) / (UMACH-DMACH) 

RETURN 

END 
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SUBPOUT I NE VPROP (T.P.D.K.H.S.U) 

COMMON /RFPR/ RF ( 1 0 ) 

.... ROUTINE TO CALCULATE THE PROPERTIES OF THE VAPOR 
.... WRITTEN 7/20/71 A MYERS 

K =1 INPUT IS T ♦ D 

K =2 INPUT IS T ♦ P 

K =3 INPUT IS T. P. ♦ 0 

IP(K.EQ.1)CALL PFNO(T.D.P) 

IF (K.E0.2) CALL OFNO (T.P.O.O) 

HOT0=RF(l) 

SOTO=RF (2) 

«FST=RF<3) 

RFHT=RF (4) 

, R =RF (5) 

AK ?RF ( 6 ) 

F1D=FIN61 (T.D) 

F10=FING1 <T»0) 

F?D=FING2 (T.D) 

F?0=F I NG2 ( T ♦ 0 ) 

SO=SOTO*CPSI (T) -CPS I (RFST ) 

HO=HOTO*CPHI (T)-CPHI ( RFHT ) 

S=SO- (R*ALOG (D*R*T ) -F 1 D*F 1 0 ) * AK 

H=HO* (T*<F1D-F10) *F2D-F20*P/D-R*T) *AK 

U=H- (P/D) *AK 

RETURN 

END 


SUBROUTINE VSND ( T * P ♦ D . K • W ) 

COMMON /RFPR/ RF ( 1 0 ) 

.... ROUTINE TO CALCULATE THE SONIC VELOCITY FOR FOLLOWING INPUT OF K 
.... WRITTEN 2/3/72 A MYERS 


K 

= 1 

INPUT 

IS 

T 

♦ P 

RFTURNS 

SONIC 

VELOCITY* W ♦ 0 


K 

= 2 

INPUT 

IS 

T 

♦ D 

RETURNS 

SONIC 

VELOCITY* W 


K 

= 3 

INPUT 

IS 

T 


RFTURNS 

W* D* 

♦ P FOR SATURATED 

VAPOR 

K 

= 4 

INPUT 

IS 

T 


RFTURNS 

W* 0* 

♦ P FOR SATURATED 

LIQUID 


AK=RF(6> 

AM=RF (7) 

IF ( (K.GT.O) .OR. (K.LT.5) ) GO TO I 
WRITE (3.300) K 

300 FORMAT ( 27H *** ERROR IN CALL VSND ***♦/♦ 

. * K MUST EQUAL 1*2* 3* OR **»/♦ 

. • K = **110) 

RETURN 

1 IF (K.EQ.2) GO TO 3 
I F ( K . GT . 2 ) GO TO 2 
CALL DFND(T.P.O.O) 

GO TO 3 

2 P= VPN ( T ) 

IF ( K • EQ • 3 ) CALL OFND (T»P*D»2) 

IF (K .E0.4) CALL DFND (T*P*D*1) 
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3 CALL CPVTD (T »D»CP»CV) 

(CP/CV)*DPOD (T.D) * ( AK*1000.00/AM) 
TE(W.LE.O. ) GO TO 4 
i rf= SORT (W) 

RETURN 

4 CONTINUE 
W = 0 

RETURN 

END 


FUNCTION CPHI (T) 

COMMON /REPR/ RE (10) 

COMMON /C I GCP/G ( 9 ) 

COMMON /SCRH/ X(40> 

.... ROUTINE TO CALCULATE INTEGRAL(CPO DT> 

.... WRITTEN 7/19/71 

R=RE(5) 

4K = RF (6) 

T?=T*T 

T3=T2*T 

T4=T3*T 

U=G(9)/T 

X<1)=-1. 00/(2. 00*T2) 

X (2)=-1.00/T 
X ( 3) = ALOG(T) 

X ( 4 ) =T 

X (5)=T2/2.00 
X (6) =T3/3.00 
X ( 7 ) = T4/4 .00 

X ( 8 ) = U* T/(FXP(U) -1.0) 

CPHI=0 
00 1 1 = 1,8 

1 CPHI =CPHI ♦ X ( I ) *G ( I ) 

C D H I =CPH I *R 

CPHI=CPHI*AK 

RETURN 

END 


FUNCTION CPIG(T) 

COMMON /RFPP/ RF (10) 

COMMON /CIGCP/G (9) 

COMMON /SCRH/ X ( 40 ) 

.... CALCULATE IDEAL-GAS HFAT CAPACITY. CP 
.... WRITTEN 7/19/71 - A MYERS 
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o=RF (5) 

AK = RF (6) 

T?=T*T 

T3=T2*T 

ll=G(9)/T 

U2=tJ*U 

X ( 1 > =1 .00/T3' 

* (2>=1.00/T2 
X (3)=1.00/T 
X (4) =1 .00 
X (5) =T 
X ( 6 ) =T2 
X (7) =T3 

X(R) s U2 * FXP (II) / (FXP <U> -1.0)««2 
CP I G = 0 .00 

do i i si, a 

1 CPIG=CPTG*X ( I ) *G ( I ) 

CPIG=CPIG»R*AK 

RETURN 

END 


FiJNCT 1 0 M CPSI (T) 

COMMON /C IGCP/G ( 9 ) 

COMMON /RFPR/ RF (10) 

COMMON /SCRH/ X ( 40 ) 

.... ROUTINE TO CALCULATE INTEGRAL (CPO/T OT) 

.... WRITTEN 7/19/71 - A MYERS 

R=RF (5) 

AK=RF (6) 

T?=T*T 
T 3=T2*T 
'J=G (9) /T 
FU = EXP (U) 

X(l)=-1. 00/(3. 00*T3) 

X (?) s-i .00/ (P.00*T2) 

X (3) s-i ,00/T 
X ( 4 ) =ALOG ( T ) 

X (5) =T 

X (6) =T2/2.00 
X (7) ST3/3.00 

X (8) su/ (EU-1 .00) -ALOG ( 1 . 00-1 .00/EU) 

CPSI=0 

DO 1 1=1*8 

1 CPSI=CPSI*X (I)*G(I) 

CPSI=CPSI*R 

cpsi=cpsi*ak 

RETURN 

END 
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FUNCTION DELC(DD.TT) 

C SET UP FOR NITROGEN 

C INPUT* TEMP=K* DFNSITY=GM/CC 

C OUTPUT* MW/CM-K 

TYPE REAL K 

DATA (AV=6.025E*23) * <K=1 .38054E-16) , (F 1 = 1 00.) . <EM=28.016) , (EOVERK= 
1118.0) * <RM=3.933E-08> » (X=1.6711) . (TC=126.26) . <DC=0.31406) 
OATA(IND=l) 

IF(IND.GT.l) GO TO 1 

C1=K**0.5/ (6. *3. 14159* (AV/EM) **0.5) 

P1=RM**2.5* ( AV/EM) **0.5*EOVERK**0.5*X 
OF 1 = 1000. /EM 
OF2=1.0 
I ND = 2 

1 »=pi*DD**0.5/TT**0.5 
FTA=VISC (OD.TT) /1000. 

RTT = ( (TT-TC) /TO **2 
«DD= ( (OD-OC) /DC) **4 
T = TT 

C CONVERTS DENSITIES TO MOLFS / LITER FOR COMPUTATION OF DERIVATIVES 

dmoles=dd*dfi 

OPT=DPOT (TT.OMOLES) *DF2*1 .01325F*6 
i)oq=DPDD (TT.DMOLES) *DF1*1,01325E*6 
0ELC02=C1*T**1 ,5*DPT**2/ (P*ETA*DD*DPD**0.5) /FI 
DFLC=DELC02 * EXP (-18.66 * BTT -4.25 * ROD) 

oflc=delc/ioo. 

RETURN 

END 


function dpdd(t.O) 
common /RFPP/'RFdO) 

COMMON /CEOS/G (4 1 ) 

COMMON /SCRH/ X (40 ) 

c 

c.... CALCULATES DP/DD of the FOUATION OF STATE 
C 

P=RF (5) 

02=D*D 

D3=D2*0 

D4=D3*D 

05=D4*D 

D6=D5*D 

07=D6*0 

D8=D7*D 

D9=D8*D 

D 1 0=D9*D 

01 1=D10*D 

o 1 2=D 1 1 *D 

013 = 012*0 

TS=SORT(T) 

t?=t*t 

T 3 = T2*T 
T4=T3*T 
Gm=G (41 ) 
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F = EXP(GM *D2) 

F1=2.00*F*GM*D 

F21=3.000*F*D2 *F 1 *D3 

F22=5.000*F*D4 *F1*D5 

F23=7.000*F*D6 *F1*D7 

F24=9.000*F*D8 *F1*D9 

F25=11.00*F*D10*F1*D11 

F26=13.00*F*D12*F1*D13 

X( 1 ) =2 • 0 0*D*T 

X( 2)=2.00*D*TS 

X( 3)=2.00*D 

X( 4 ) =2 • 0 0*D/T 

X( 5 ) =2 • 0 0*D/T2 

X( 6)=3.00*D2*T 

X( 7)=3. 00*02 

X( 8 ) =3 . 00*D2/T 

X( 9) =3.00*D2/T2 

X(10)=4.00*D3*T 

X ( 1 1 ) =4 , 0 0*D3 

X ( 12) =4.00*D3/T 

X ( 13) =5.00*04 

X ( 14) =6. 00*D5/T 

X ( 15) =6.00*D5/T? 

X ( 16) =7. 00*D6/T 
X ( 17) =8.00*D7/T 
X ( 18) =8. 00*D7/T2 
X ( 19) =9.00*D8/T2 
X (20)=F21/T2 
X (21 ) =F21/T3 
X (22) =F22/T2 
x (23) =F22/T4 
X (24>=F23/T2 
x (25) =F23/T3 
X ( 26 ) =F24/T2 
X (2 7)=F24/T4 
X (28) =F25/T2 
x (28) =F25/T3 
x ( 30 ) =F26/T2 
X ( 3 1 ) =F26/T3 
X (32) =F26/T4 
'J = 32 
0PDD = 0 
00 1 K=1*N 

1 OPDO = DPDO*X (K) *G(K), 
OPD0 = DPD04.R*T 
RETURN 

EmO 


FUNCTION DPDT ( T * 0 ) 
COMMON /CEOS/G (41) 
COMMON /RFPR/ PF (10) 
COMMON /SCPH/ X ( 4 0 ) 
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C 

c.... CALCULATES dp/dt of the equation of state 
c 

R=RF (5) 

02=0*0 

03-02*0 

04=03*0 

05=04*0 

06=05*0 

07=06*0 

08=D7*D 

09=08*0 

0 1 0=D9*O 

0 1 1=0 1 0*D 

0 1 2=D 1 1 *0 

013=012*0 

TS=SQRT(T) 

T?=T*T 
T 3 = T 2 * T 
T4=T3*T 
T5=T4*T 
r;M=G (41 ) 

F = EX p ( GM *02) 

X( 1 ) =0? 

X( 2) =02/ ( 2 • 00*TS ) 

X( 3 ) =0 

X( 4 ) =-02/T2 

X( 5)=-?*00*D2/T3 

*( 6) =03 

X< 7 ) = 0 

X( 8 ) =-D3/T? 

X( 9) =-2.00*D3/T3 
X (10) =04 
X (11)=0 
X (12)=-D4/T2 
x ( 1 3 ) =0 
x (14) =-D6/T2 
X ( 15) =-2‘.00*D6/T3 
X ( 16) =-07/T2 
X (17)=-08/T2 
X (18)=-2 # 00*08/T3 
X ( 19) =-2.00*DQ/T3 
X (20)=-2#00*D3*F/T3 
X(21)=-3.00*D3*F/T4 
X(22)=-2.00*D5*F/T3 
X (23)=-4.00*D5*F/T5 
x (24)=-2.00*D7*F/T3 
X (25)=-3.00*D7*F/T4 
X (26)=-2.00*D9*F/T3 
X ( 27 ) =*4 . 0 0*D9*F/T5 
X (28)=-2.00*011*F/T3 
X(29)=-3,00*D11*F/T4 
X ( 30 ) =-2.00*0 1 3*F/T3 
x (31) =-3*00*013 *F/T4 
X (32)=-4.00*D13*F/T5 
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N=32 
DPDT=0 
DO 1 K=ltN 

1 DPDT=DPDT*X(K)*G(K) 
DPDT=DPDT ♦R*D 
RFTURN 
END 


FUNCTION DSATL(T) 

COMMON /CRPR/CR ( 3 ) /CSL/G<7> 

COMMON /SCRH/ B ( 40 ) 

.... THIS FUNCTION SUPPLIES AN APPROXIMATE VALUE FOR THE 
.... density OF THF saturated LIQUID 

.... WRITTEN 9/20/71 A Myf'RS 

TC=CR (3) 

X= (TC-T) /TC 

x?=x*x 

X3=X*X2 

X4=X*X3 

X5=X*X4 

B ( 1 ) =1 .00 

H < 2 ) =X 

B (3) =X2 

B(4)=X3 

B(5)=X4 

B (6) =X5 

B ( 7 ) =ALOG ( X ) 

DSL = 0 

DO l 1=1,7 
1 OSL=DSL*B ( 1 1 *G ( I ) 

OSATL=DSL 

RETURN 

END 


FUNCTION DSATV(T) 

COMMON /CRPR/CR ( 3 ) /CS V/G ( 7 ) 

COMMON /SCRH/ B(40) 

C 

C.... THIS FUNCTION SUPPLIES AN APPROXIMATE VALUE FOR THE 
C.... DENSITY OF THF SATURATED VAPOR 
C 

C.... WRITTEN 9/20/71 A MYERS 
C 

TC=CR<3) 
x= (TC-T) /TC 
X2=X*X 
X3=X*X2 
X4 =x*X3 
X S = X * X4 
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B ( 1 ) = 1 . 00 

B(2)=X 

B ( 3 ) =X2 

B(4)=X3 

B(5)=X4 

B(6)=X5 

B ( 7) =ALOG ( X ) 

DSV=0 

DO 1 1*1.7 
1 OSV=DSV+BCI)*G(I> 
DSV = EXP(DSV) 
DSATV=DSV 
RETURN 
END 


FUNCTION DSFND(SOT.TT.DD) 

T = TT 
0=0D 

00 10 1 * 1.100 

CALL PROP (T .P.0.2.H.S.U) 

IF (A8S(S-SOT)-1.F-07*SOT) 11.11.1 
1 OUP=D* » 00 1 *D 
DDN=D-.001*D 

CALL PROP (T .P.DUP.2.H.SUP.U) 

CALL PROP ( T * P.DDN. 2.H , SDN.U) 
0=D-(S-S0T)/( <SUP-SDN)/(OUP-DDN) ) 

10 CONTINUE 
0=0 

11 OSFN0=0 
RETURN 
END 


FUNCTION FINGl(T.O) 

COMMON/CEOS/G (41 ) 

COMMON /SCRH/ X (40 ) 

.... ROUTINE TO CALCULATE INTEGRAL ( (R/D-l/D**2 (DP/DT) ) OD) 

.... WRITTEN 7/1B/71 A MYERS 

D2=D*D 

03=02*0 

04=03*0 

05=D4*D 

06=05*0 

D7 = D6*0 

08=07*0 

09=08*0 

010=09*0 
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TS=SQRT(T> 

T2=T*T 
T3=T2*T 
T4=T3*T 
T5=T4*T 
GM=G(41 ) 

F = EXP (GM *02 ) 

S1=F/(2.00*GM) 

G2= (F*02-2. 00*61 ) /(2.00*GM) 
63= (F*D4-4.00*G2) /(2.00*GM) 
64= (F*D6-6.00*G3>/ (2.00*GM) 

65 = (F*D8-8. 00*64 ) / (2.0 0*GM) 
66* (F*010-10.00*G5) / (2.00*GM) 
X( l)=-0 

X( 2>=-D/(2.00*TS) 

X( 3) =0 
X ( 4 ) =*D/T2 
X( 5>=2.00*D/T3 
X( 6 ) =-D2/2 .00 
X( 7) =0 

X( 8 ) =62/ ( 2 • 0 0*T2 ) 

X( 9>=02/T3 
X ( 1 0 ) =-D3/3. 00 
X ( 1 1 ) =0. 

X (12)=D3/(3.00*T?) 

X ( 1 3 ) =0 

X ( 14) =05/(5. 00*T2) 

X ( 1 5 ) * 2.00*05/ (5. 00*T3) 

X ( 16) *06/ (6.00*T2) 

X (17) =07/ (7.00*T2> 

X ( 18)=2.0 0*D7/ (7.0 0*T3> 

X (19)=08/(4.00*T3) 

X (20) =2.00*G1/T3 
X (21 ) =3.00*G1/T4 
X (22) =2.00*G2/T3 
X (23) =4.00*62/75 
X (24)=2.00*G3/73 
X (25)=3. 00*63/74 
x (26) =2.00*64/T3 
X (27) =4.00*G4/T5 
X (28) =2.00*65/T3 
X (29) =3.00*65/T4 
X (30) =2.00*66/T3 
X <31)=3.00*G6/T4 
X ( 32 ) =4 . 0 0*66/T5 
FING1=0 
00 1 1=1.32 

1 FING1=FING1>G(I)*X(I) 

RETURN 

END 


FUNCTION FING2 (T.O) 
COMMON/CEOS/G (41) 
COKMON /SCRH/ X(40) 



o o o o o 


APPENDIX C 


.... ROUTINE TO CALCULATE INTEGRAL ( (P/D**2-RT/D) DO) 

.... WRITTEN 7/18/71 - A MYERS 

02=D*0 

03=D2*D 

04=D3*0 

05=04*0 

06=D5*0 

07=06*0 

08=07*0 

09=08*0 

010=09*0 

TS=SQRT (T) 

T?=T*T 

T3=T2*T 

T4=T3*T 

T5=T4*T 

GM=6(41) 

F = EXP(GM *D2 ) 

G1=F/(2.00*GM) 

PRs(F*02-2.00*G1 ) / <2.00*GM> 

G3= (F*D4-4.00*G2) /(2.00*GM) 

G4= (F*06-6. 00*63) / (2.00*6M) 
G=i=<F*OB-8.00*G4>/(2.00*GM) . 

G6=(F*D10-10 .00*65 )/(2.00*GM) 

X( 1 > =0*T 
X( 2 ) =0*TS 
X( 3) =0 
X( 4 ) =D/T 
X ( 5) =0/T2 
X( 6)=02*T/2.00 
X( 7 ) =02/2 .00 
X ( 8 ) =02/ ( 2 • 00*T ) 

X ( 9)=D2/(2.00*T2) 

X(10)=03*T/3.00 
X ( 1 1 ) =03/3 , 00 
X(12)=D3/(3.00*T) 

X ( 1 3) =04/4 .00 
X (14)=05/(5.00*T) 

X(15)=05/(5.00*T2> 
x (16)=06/(6.00*T) 

X ( 17 ) =07 / ( 7. 00*T ) 

X (18>=D7/(7.00*T2) 

X ( 19) =08/ (8.00*T2) 

X (20) =G1/T2 
X (21)=G1/T3 
X ( 22 ) =G2/T2 
X (23) =G2/T 4 
X (?4)=G3/T2 
X (25)=G3/T3 
x (26)=G4/T2 
X ( 27 ) =G4/T4 
X ( 28 ) =G5/T2 
X ( 29 ) =G5/T 3 
X ( 30 ) =G6/T2 
x (31) =G6/T3 
X ( 32 ) =66/T4 
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APPENDIX C 


FTNG2=0 
00 1 1=1,32 

1 FING2=FING2*G(I)*X (I) 
RFTURN 
FND 


FUNCTION F I NG3 ( T , 0 ) 

C0MM0N/CE0S/G(41 ) 

COMMON /SCRH/ X(40) 

C 

C.... ROUTINE TO CALCULATE I NTF GRAL < ( T/D**2 ) * ( D2P/DT2 > DO) 
C 

C.... WRITTEN 7/16/71 - A MYERS 
C 

02=0*0 

D3=D2*D 

04=03*0 

D5=D4*0 

06=D5*D 

07=06*0 

08=D7*D 

09=08*0 

010=09*0 

TF=SORT (T) 

T?=T*T 
T 3 = T 2 * T 
T4=T 3*T 
T5=T4*T 
GM=G(41) 

F = EXP ( GM *D2 ) 

0 1 =F / ( 2 • 0 0*GM ) 

G? = (F*02-2.00*G1 ) / (2.00*GM) 

G3= (F*04-4. 00*G2) / (2. 00*GM) 

G4= (F*D6-6.00*G3) / ( 2 • 0 0*GM ) 

G5= (F*D8-8.00*G4) / (2.00*GM) 

G8= (F*D 10-10. 00*G5) / <2.00 *GM) 

X ( 1 ) =0 

X< ?)=-D/(4,00*TS) 

X ( 3 ) =0 

X( 4) =2.00*D/T2 
X( 5) =6. 00*0/T3 
X( 6 ) = 0 
X ( 7) =0 
X( 8 ) =02/T2 
X( 9 ) =3 . 0 0*O2/T3 
X(10)=0 
X(11)=0 

X ( 12) = (2.00*03) / <3.00*T2) 

X ( 1 3 > =0 

X ( 14) = (2.00*D5) / (5.00*T2) 

X ( 15) = (6.00*05) / (5.00*T3) 
x ( 16) =06/ (3.00*T2) 

X (17) = (2.0 0*D7)/(7.00*T2) 

X (18) = (6.00*07)/ ( 7 * 0 0*T3 ) 

X ( 19) = (3.00*08) / (4.0 0*T3) 

X (20) =6.000*G1/T3 
X (21 ) =12.00*G1/T4 '• 
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X (22)=6.000*G2/T3 
X (23)=20.00*G2/T5 
X(24)=6.000*G3/T3 
X(25)=12.00*G3/T4 
x (26>=6,000*G4/T3 
x(27)=20.00*G4/T5 
X (28>=6.000*G5/T3 
X(29)=12.00*G5/T4 
X (30) =6.000*G6/T3 
X(31)=12.00*G6/T4 
X(32)=20.00*G6/T5 
FING3 = 0 
On 1 1=1,32 

1 PINIG3=FING3*G(I)*X(I) 
RETURN 

F i'JD 


C 

c 

c 

c 

c 

c 

1 


FUNCTION VISC ( DO * TT ) 

OTMENSION 7(18). X<11> 

DATA ( Z = -6.8939127475E*01 


3.52261 18983E*00 


-6.8357539823F-02 
3.6093309138F-09 
-1.0717599406E-19 
3.853077101 1F-03 
8.905971 1315E-10 
-2.3084044942E-20 


1 .5832717315E-03 . 

-2.5555598476E-12 . 

7.4165322904E+01 , 

8.0133713668E-04 , 

-5.3779372664E-13 ■ 


-2.6418423047E-06 
8.5635041641E-16 
-1 .5834400475E+00 
-8.9203123846E-07 
1.7398277309E-16 


DATA (X = 2.3083514362E-1 » -9. 36362071 7 IE-1 ,9. 0339 186452E + 0, -4. 1832067 
163E-H .1 ,0R97627893E*2.-1 . 29 1 3856376E ♦ 2 , 5 . 978204991 3E* 1 ) 

J = 9 $ K = 2 

GO TO 100 


ENTRY THERM 
J = 0 * K = 1 

100 CONTINUE 

VISC = 0 S T = TT $ D = DD 
00 200 1=1,9 

VISC = VISC ♦ T**(I-3)*Z(I*J) 

200 CONTINUE 

GO TO (300,400) ,K 
300 CONTINUE 

VISC = VISC + .21959022190 * D ♦ . 063873706990* ( E XP ( 3 . 6*0 ) - 1. ) 

RETURN 

400 CONTINUE 

IF (D.GT.0.8) GO TO 501 
DO 500 1=1,7 

VISC * VISC ♦ X(I) * D**I 
500 CONTINUE 
RETURN 

501 VISC=VISC-2.5781990818E*3+9.4784808659E*3*D-1 . 1622926973E*4*D**2* 
14.756948538E*3*0**3 
C INPUT 0=GM/CC, T=DEG K 

C OUTPUT ETA=MILI GM/CM-SEC L AMBO A=MW/CM-K 

RETURN 
END 


M2 


o o o o o 


APPENDIX C 


FUNCTION! VPN(T) 

COMMON/CVPN/G( 11 ) /CRPP/CP ( 3 ) 
COMMON /SCRH/ X(40) 

.... CALCULATE THE VAPOR PRESSURE 

.... WRITTEN 7/22/71 A MYFPS 

TC=CR ( 3 ) 

A = G ( 1 1 ) 

T2=T*T 
T3=T*T2 
T 4 = T *T3 
TS = T*T 4 
T G = T *T5 
X ( 1 ) =1 .00/T 
X (2) =1 .00 
X (3) =T 
X (4) =T2 
X (5) =T3 
X (6)=T4 
X (7) =T5 
X (8) =T6 
x (9) =ALOG (T) 

X ( 10) = (TC-T) **A 

d=0 

DO 1 1=1.10 

1 P = P + X ( I ) *G ( I ) 

P=EXP (P) 

VPN = P 

RETURN 

END 
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Sample Output 


ENERGY AND POWER REQUIREMENTS FOR FAN-DRIVEN CRYOGENIC WIND-TUNNEL 
ISEN TROPIC COMPRESSION IN NITROGEN 


TEST-SECTION MACH NUMBER 
FAN PRESSURE RATIO 
FAN PRESSURE# OUTLET 
FAN PRESSURE# INLET 


1.20 
1.2000 
5.0000 ATM 
A. 1667 ATM 


E - FAN ENERGY 
P • FAN POWER 




REAL 

GAS 


— 

IDEAL 

GAS 


REAI 

L/ IDEAL VAI 

LUES 

TT# IN 

T T# OUT 

MDOT 

E/MASS 

P/AREA 

TT# OUT 

MDOT 

E/MASS 

P/AREA 

MDOT 

E /HASS 

P/AREA 

K 

K 

KGM/S-M2 

J/KGM 

J / S— M2 

K 

KGM/S-M2 

J/KGM 

J/S-M2 




294.2 

310.0 

1145.5 

16332.7 

. 1871E+08 

309.9 

1143.8 

16341.6 

.1869E+08 

1.0015 

.9995 

1.0010 

284.7 

300.0 

1164.7 

15799.5 

• 1840E +08 

299.9 

1162.7 

15814.1 

•1839E+08 

1.0017 

.9991 

1.0008 

275.2 

290.0 

1184.9 

15265.9 

•1809E+06 

289.9 

1182.6 

15286.7 

•1808E+08 

1.0020 

.9986 

1.0006 

265.7 

280.0 

1206.2 

14731.7 

• 1777E + 08 

279.9 

1203.5 

14759.2 

•1776E+08 

1.0022 

.9981 

1.0004 

256.2 

270.0 

1228.7 

14197.0 

.1744E+08 

269.9 

1225.6 

14231.8 

« 1744E+08 

1.0025 

.9976 

1.0001 

246.7 

260.0 

1252.6 

13661.5 

•1711E+08 

259.9 

1249.0 

13704.4 

.1712E+08 

1.0029 

.9969 

.9997 

237.2 

250.0 

1277.9 

13125.3 

•1677E+08 

249.9 

1273.7 

13176.9 

.1678E+08 

1.0033 

.9961 

.9993 

227.7 

240.0 

1304.9 

12588.2 

• 1643E+0 8 

239.9 

1300.0 

12649.5 

0 1644E +08 

1.0037 

.9952 

.9989 

218.2 

230.0 

1333.7 

12050.1 

• 1607E +0 8 

229.9 

1328.0 

12122.1 

. 1610E+08 

1.0043 

.9941 

.9983 

208.8 

220 .0 

1364.5 

11510.7 

. 1571E+08 

219.9 

1357.9 

11594.6 

•1574E+08 

1.0049 

.9928 

.9976 

199.3 

210.0 

1397.7 

10969.9 

• 15 3 3E +0 8 

209.9 

1389.8 

11067.1 

•1538E+08 

1.0056 

.9912 

.9968 

189.8 

200.0 

1433.5. 

10427.5 

. 149 5E + 08 

199.9 

1424.2 

10539.6 

.1501E+08 

1.0065 

.9894 

.9958 

180.3 

190.0 

1472.3 

9883.0 

. 1455E+08 

189.9 

1461.2 

10012.1 

•1463E+08 

1.0076 

.9871 

.9946 

170.8 

180.0 

1514.7 

9336.0 

• 1 41 4E +0 8 

179.9 

1501.3 

9484.6 

•1424E+08 

1.0089 

.9843 

.9931 

161.3 

170.0 

1561.1 

6786.1 

. 1372E+08 

169.9 

1544.9 

8957.0 

•1384E+08 

1.0105 

.9009 

.9912 

151.8 

160.0 

1612.5 

823Z.3 

•1327E+08 

159.9 

1592.5 

6429.5 

.1342E+08 

1.0125 

.9766 

.9889 

142.3 

150.0 

1669.7 

7673.7 

. 12 8 IE + 08 

149.9 

1644.8 

7901.9 

. 1300E+08 

1.0151 

.9711 

.9858 

132.8 

140.0 

1734.2 

7108.6 

•1233E+Q8 

139.9 

1702.7 

7374.2 

• 1256E+08 

1.0185 

.9640 

' .9818 

123.3 

130.0 

1807.8 

6534.8 

.1181E+08 

129.9 

1767.0 

6846.6 

.1210E+08 

1.0231 

.9545 

.9765 

113.8 

120.0 

1893.5 

5948.5 

.1126E+08 

119.8 

1839.3 

6318.9 

.1162E+08 

1.0294 

.9414 

.9691 

109.0 

115.0 

1942.2 

5648.9 

• 1097E+08 

114.8 

1879.0 

6055.1 

. 1138E+08 

1.0336 

.9329 

.9643 

104.3 

110.0 

1995.8 

5343.4 

• 1066E+08 

109.8 

1921.3 

5791.2 

•1113E+08 

1.0388 

.9227 

.9584 

99.5 

105.0 

2055.5 

5030.5 

. 1034E+08 

104.8 

1966.6 

5527.4 

• 10 B7E+08 

1.0452 

.9101 

.9512 


SATURATION CONDITIONS PEACHED AT THROAT MACH NUMBER 


Output variables: TT,IN, stagnation temperature at fan inlet, K; TT,OUT, stagnation temperature at 

fan outlet, K; MDOT, mass-flow rate per unit throat area, kgm/sec-m 2 : E/MASS, fan energy per unit 
mass of flow, J/kgm; P/AREA, fan power per unit throat area, J/see-m^. 
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Figure 2.- Analytical model of tunnel. 



PROGRAM INPUTS 



Figure 3.- Flow chart for isentropic power requirements program. 
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(a) Pt,2 = 1 -° atm * 

Figure 5.- Energy for isentropic compressions of nitrogen (relative to ideal-gas values) 
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Figure 7.- Variation of isentropic compression energy for nitrogen with stagnation pressure 

(relative to ideal-gas values), r = 1.025. 
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Figure 8.- Concluded 
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(a) p. 5 '= 1.0 atm. 

Figure 9.- Relative isentropic power values for various stagnation temperatures and pressures. 
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(b) Pt t 2 = 3-0 atm. 


Figure 9.- Continued. 
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Figure 11.- Estimates of energy for isentropic compressions of nitrogen compared with 

exact values. pt f 2 = atm » r = 1*20. 
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Figure 12.- Approximation for real-gas mass-flow rates of nitrogen (isentropic flow). 
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Figure 15.- Flow chart for real-gas calculations of fan parameters. 
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